arXiv:hep-lat/0512017v2 12 Apr 2006 


DESY 05-250 
MS-TP-05-33 
SFB/CPP-05-83 


Numerical simulations with two flavours of 
twisted-mass Wilson quarks and DBW2 gauge action 

F. Farchioni“, P. Hofmann®, K. Jansen^, L Montvay^, 

G. Miinster®, E.E. Scholz^* L. Scorzato^, A. Shindler^, 

N. Ukitad C. UrbacB’®| U. Wenger^§ L Wetzorke^ 

“ Universitat Munster, Institut fiir Theoretische Physik, 
Wilhelm-Klemm-Strasse 9, D-48149 Munster, Germany 
^ NIC/DESY Zeuthen, Platanenallee 6, D-15738 Zeuthen, Germany 
Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, D-22603 Hamburg, Germany 
^ Institut fiir Physik, Humboldt Universitat zu Berlin, D-12489 Berlin, Germany 
® Freie Universitat Berlin, Institut fiir Theoretische Physik, 

Arnimallee 14, D-14196 Berlin, Germany 


Abstract 

Discretisation errors in two-flavour lattice QGD with Wilson-quarks and DBW2 
gauge action are investigated by comparing numerical simulation data at two values of 
the bare gauge coupling. Both non-zero and zero twisted mass values are considered. 

The results, including also data from simulations using the Wilson plaquette gauge 
action, are compared to next-to-leading order chiral perturbation theory formulas. 

*Present address: Physics Department, Brookhaven National Laboratory, Upton, NY 11973 USA 
^Present address: Institute of Physics, University of Tsukuba, Tsukuba, Ibaraki 305-8571, Japan and 
Department of Physics, The University of Tokyo, Kongo 7-3-1, Bunkyo-ku, Tokyo 113-0033, Japan 

^Present address: Theoretical Physics Division, Dept, of Mathematical Sciences, University of Liverpool, 
Liverpool L69 3BX, UK 

^Present address: Institute for Theoretical Physics, ETH Ziirich, CH-8093 Zurich, Switzerland 


1 



1 Introduction 


The singular point of QCD at vanishing quark masses is distorted in Wilson-type lattice 
formulations: as a result of lattice artefacts, in the region of small quark masses an 
extended phase structure is developed. This phase structure can be predicted and 
classified in Chiral Perturbation Theory (ChPT) ^ if lattice artefacts are taken into 
account j^j- If, in addition to the usual quark mass parameter, a twisted quark mass 
is introduced 011 ] then in the plane of untwisted and twisted quark mass a first order 
phase transition line with second order endpoints appears. Depending on the sign of 
the leading term representing lattice artefacts, the first order phase transition line is 
either on the untwisted quark mass axis (“Aoki-phase scenario” 0) or perpendicular 
to it (“normal scenario”) 000 . 

In numerical simulations it pays off to try to reduce lattice artefacts at fixed (non¬ 
vanishing) lattice spacing by an appropriate choice of the lattice action. An important 
issue in this respect is to bring the phase structure at small quark masses as close 
as possible to the point-like singularity appearing in the continuum limit. In fact, 
the strong first order phase transition observed earlier in numerical simulations with 
Wilson-type quarks 00[in] presents a serious obstacle for QCD simulations with light 
quarks. 

In previous work we systematically investigated the phase structure of lattice QCD 
with twisted-mass Wilson-type quarks (for a recent review see [1 1 1 ). In Ref. m we have 
shown that at lattice spacings near a ~ 0.2 fm the phase structure with Wilson-quarks 
and Wilson-plaquette gauge action is consistent with the “normal scenario” of ChPT. 
This differs from the situation in the strong coupling regime, where the “Aoki-phase 
scenario” has been previously observed HSj. 

A consequence of the “normal scenario” is that for fixed gauge coupling {(3) the 
mass of charged pions have a positive lower bound (m™"). The numerical simulation 
data in Ref. m have shown that this lower bound is at a ~ 0.2 fm quite high, namely 
about 600 MeV. Such a high lower bound would prohibit the study of light quarks. 
Therefore, an important question is the behaviour of this lower bound as a function of 
the gauge coupling (or lattice spacing) towards the continuum limit. In a subsequent 
paper it has been shown m that, as expected, the lower bound becomes clearly smaller 
for decreasing lattice spacing. Its decrease in the range 0.20 fm > a > 0.14fm is roughly 
consistent with the prediction of next-to-leading-order (NLO) ChPT 0000II3, 
namely oc a (at a/r = 0). A minimal pion mass of ~ 300 MeV is estimated 
to occur near a ~ 0.07 — 0.10 fm, but this estimate is rather uncertain and has to be 
checked in future simulations if the Wilson gauge action ought to be used. The question 
arises whether one could lower m™*"' by a suitable change of the lattice action. 

An early observation by the JLQCD Collaboration has been 0 that the strength 
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of the first order phase transition near zero quark mass is sensitive to a change of 
the gauge action. Following this hint, we have shown in a previous paper m that 
combining two flavours {Nf = 2) of Wilson-quarks with the DBW2 gauge action |17j 
leads to a phase structure near zero quark mass with substantially weaker first order 
phase transition. As a consequence, the minimal pion mass is at least by a factor of 
two lower compared to the plaquette gauge action at similar lattice spacings. 

This implies that numerical simulations with light quarks become possible on 
coarser lattices and hence with much less computational costs if the DBW2 gauge 
action is used. Of course, for the choice of the gauge action also other criteria may be 
relevant. For instance, it has been reported in quenched studies HHHH that in some 
quantities strong scale breaking effects appear if the DBW2 action is used. Another 
problem could be the late convergence of lattice perturbation theory, implied by the 
results of the QCDSF Collaboration |2nj . 

In general, the question of the scaling behaviour of the results obtained by a given 
lattice action is very important. In case of the Wilson twisted-mass formulation of 
lattice QCD it has been shown m that the leading lattice artefacts are of O(a^) if 
the bare quark masses are appropriately tuned. Detailed investigations have shown 
[H usi nn that in the quenched approximation excellent scaling behaviour can be 
achieved, indeed, also at light quark masses. The same question in the full theory with 
dynamical quarks is obviously very important. 

In the present paper we perform first exploratory scaling tests for the combination 
of Wilson-fermion lattice action with the DBW2 gauge action by comparing numerical 
simulation data at two values of the gauge coupling, namely (3 = 0.67 and (3 = 0.74. 
We consider data points with both vanishing and non-vanishing value of the twisted 
mass. Moreover, since one can extract useful information on multiplicative renormali¬ 
sation factors from the dependence of matrix elements on the twist angle in the plane 
of untwisted and twisted quark mass, we exploit this method and derive from our 
simulation data the values of Zy, Za and ZpjZs- In addition, we compare the NLO- 
ChPT formulas of Refs, [hi l6l 171 ITbl to the results of the numerical simulations. For 
comparison, ChPT fits of the data obtained by the Wilson plaquette gauge action na 
are also considered. 

The outline of the paper is as follows: in the next section, after specifying the lattice 
action and the simulation algorithms, the numerical simulation runs are discussed and 
some scaling tests are presented. Section |21 is devoted to a detailed description of the 
results on the twist angle in the plane of untwisted and twisted quark mass together 
with an explanation how the aforementioned multiplicative renormalisation Z-factors 
can be determined. The knowledge of the twist angle and Z-factors makes it possible 
to obtain results on physical quantities, such as the quark mass and the pion decay 
constant. In Section^the ChPT fits of the data with DBW2 gauge action are presented. 
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Section 13 contains a discussion and a summary. In an Appendix alternative chiral fits 
of the DBW2 data are shown and compared to similar ChPT fits of Wilson plaquette 
data. 


2 Numerical simulations 

The lattice action and simulation algorithms are defined here for the reader’s conve¬ 
nience. The notations are similar to those in Ref. [El. 


2.1 Lattice action and simulation algorithms 

We apply for quarks the lattice action of Wilson fermions, which can be written as 


Sq = '^{{XxW + ilbT:ia^i]xx) (xx+iiUxp.[r + ly\Xx) 

Here the (“untwisted”) bare quark mass in lattice units is denoted by 

, 1 

/is = 0772-0 + 4r = — , 

2k 


( 1 ) 


( 2 ) 


r is the Wilson-parameter, set in our simulations to r = 1, amo is another convention for 
the bare quark mass in lattice units and k is the conventional hopping parameter. The 
twisted mass in lattice units is denoted here by o/x. (This differs from the notation in 
d where /x has been dehned without the lattice spacing factor a in front.) Ux^ G SU(3) 
is the gauge link variable and we also defined Ux-jj. = ^ and 7 _^ = —7/^. 

For the SU(3) Yang-Mills gauge field we apply the DBW2 lattice action which 
belongs to a one-parameter family of actions obtained by renormalisation group con¬ 
siderations. Those actions also include, besides the usual (1x1) Wilson loop plaquette 
term, planar rectangular (1 x 2) Wilson loops: 





l--Ret/i;J 


+ ci ^ 


l--ReC/i;J 


( 3 ) 


with the condition cq = 1 — 8ci. For the DBW2 action we have ci = —1.4088. 

For preparing the sequences of gauge configurations two different updating algo¬ 
rithms were used: the Hybrid Monte Carlo (HMC) algorithm [E] with multiple time 
scale integration and mass preconditioning as described in EH and the two-step multi¬ 
boson (TSMB) algorithm m which has been tuned for QCD applications following 

Eniini. 
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2.2 Simulation parameters and a first scaling test 

In our numerical simulations we considered two values of the gauge coupling, namely 
{3 = 0.67 and (5 = 0.74. The simulations at the lower /3-value have been performed 
on a 12^ . 24 lattice as in cni- The higher /3-value (/3 = 0.74) was chosen in such a 
way that the physical volume of the 16^ • 32 lattice remains approximately the same, 
that is a(/3 = 0.74) ~ | a(/3 = 0.67). The value of the lattice spacing was defined by 
extrapolating the Sommer scale parameter in lattice units ro/a 1301 to zero quark mass 
and assuming tq = 0.5 fm. The simulation parameters and the amount of statistics are 
specihed in Tabled 

As Tabled shows, both zero and non-zero twisted mass points were simulated. The 
non-zero values of the twisted mass were also chosen according to the assumed scale 
ratio, that is o/x(/3 = 0.74) = | o/x(/3 = 0.67) = 0.0075. In other words, the bare 
twisted mass /r is kept (approximately) constant. 

In several points of the parameter space simulation runs have been performed with 
both the HMC and the TSMB updating algorithms. Having run the two algorithms in 
the same points allowed to compare their performance. It turned out that the optimised 
HMC algorithm of Ref. m is substantially faster than TSMB. For instance, in long 
runs at the simulation point (A) (16^ • 32 lattice, /3 = 0.74, k = 0.1580, o/r = 0) HMC 
with multiple time scale integration and mass preconditioning is almost by a factor 
of 10 faster. Therefore, in the majority of simulation points the hnal data analysis is 
based on HMC runs. Results from TSMB updating were only used in the runs of the 
first part of Tabled (those at /3 = 0.67 and a/r = 0). Even if results with both updating 
algorithms were available in several other points, in the final analysis we never mixed 
results from different updating procedures. 

The results for some basic quantities are collected in Tables |21 and 01 The pseu¬ 
doscalar meson (“pion”) mass am-j^ is obtained from the correlator of the charged 
pseudoscalar density 

Pt = Xx^lbXx (4) 

where r± = ri ± iT 2 . In case of the vector meson (“p-meson”) mass amp, for generic 
values of the bare untwisted and twisted quark mass, the correlators of both vector 
iy^p) and axialvector (A“^) bilinears of the x-fields can be used: 

= Xx^Ta'JpXx , = Xx\TallilbXx (« = 1,2). (5) 

The reason is that the physical vector current is, in general, a linear combination of 
Vy and A%p (see Section El). In a given simulation point we determined amp from the 
correlator possessing the better signal. 

In Table El the values of the bare (untwisted) PC AC quark mass am^^'^ are also 
given. It is defined by the PCAC-relation containing the axialvector current in 
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m and the pseudoscalar density P^\ 


am 


PCAC 

X 


{d;A+^p-) 


( 6 ) 


Here d* denotes, as usual, the backward lattice derivative. 

Besides , Table |31 also contains the values of the bare “untwisted” pseu¬ 

doscalar decay constant u/^tt defined by 


af^.^ = (am^) ^(0|H+^oo|7r ). 


( 7 ) 


The relation of the bare (untwisted) quantities and to the corresponding 

physical quantities will be discussed in the following section. 

The squared ratio of the pion mass to the p-meson mass is plotted in Figure ^ as 
a function of both of which are expected to be approximately proportional 

to the quark mass for small quark masses. (This holds if the effect of the “chiral 
logarithms” is negligible in the quark mass depedence of and if tq is approximately 
constant near zero as a function of the quark mass.) The straight line in the figure 
connects the origin and the point with the physical values m^^ = 140 MeV, rup = 
770 MeV and vq = 0.5 fm. As the figure shows, in this plot there are observable 
scale breaking effects between f3 = 0.67 and f3 = 0.74, but the (5 = 0.74 points are 
already close to the continuum expectation. Within the (large) statistical errors there 
is no noticeable difference between the points with vanishing and non-vanishing twisted 
mass. (According to Table El the twisted mass values are given by = 0.02845(68) 
and ro/i = 0.0283(15) for (5 = 0.67 and (3 = 0.74, respectively.) 


3 Twist angle and renormalisation factors 

3.1 Twist angle 

In this section we discuss the determination of the twist angle oj. For given afi) this 
is defined as the rotation angle relating twisted-mass QCD (TMQCD) to the physical 
theory QCD. An important point is that the connection can be made only after (lattice) 
renormalisation of the theory. The renormalisation of the local bilinears in the Wilson 
twisted-mass formulation is therefore involved. Some of the arguments of this section 
were already discussed in previous publications of this collaboration [Sunni- 

Following |32j we operationally define isnnni the twist angle to as the chiral rota¬ 
tion angle between the renormalised (physical) chiral currents and the corresponding 
bilinears of the twisted formulation. We denote with and the physical vector 
and axialvector currents, while and are the bilinears of the y-fields defined in 
Eq. ([Si). In order to establish the correspondence with the physical currents, the bilin¬ 
ears of the x-fields have to be properly renormalised. This is obtained, as in QCD, by 
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multiplying them by the respective renormalisation constants Zy and Z/^. In a mass 
independent scheme these are functions of /3 alone and coincide with the analogous 
quantities in Wilson lattice QCD for the same value of (3. So the relation reads: 

^ 4 . = ZyV^^ cosIV +eabZAAl^s\nu: , (8) 

= ZaA’^^^cosuj +eahZvV^^s\niv (9) 

where only charged currents are considered (a=l, 2) and €ab is the antisymmetric unit 
tensor. 

The conserved vector current of the y-fields 

\ {xx+f^TaUxt^ilf, + r)xx + Xx^aU^^ijf, - r)Xx+ti) ( 10 ) 

satisfies by construction the correct Ward-Takahashi identity of the continuum. In this 
case the Formulas (©, ® apply with Zy replaced by 1, in particular 

Ax^ = ZaAx^ cosiv + Cab 144 sinw . (11) 

In practical applications it is useful to define two further angles ivy and uja'- 

ivy = arctan(ZAZY^ taniv) , oja = arctan(ZyZ^^ taniv) . (12) 


In terms of cvy, cva Eqs. Q and © read 

144 = JVv (cos (VyV^^^ + Cab sin (Vy A^^^) , (13) 

Axf, = AfA (cos (vaA;^^ + Cab sin (vaVx^). (14) 


The unknown multiplicative renormalisations are now contained in an overall factor 
(X = V,A): 


Mx = 


Zx 


cos cJxVl + tan uja 

From the definition m it follows 


(15) 


( 16 ) 

- 4 ^ = \ /tan ivy/ tan IVA ■ ( 17 ) 

Zy ^ 

As already proposed in isDuni, we determine the twist angle iv by imposing parity- 
restoration (up to 0(a) precision) for matrix elements of the physical currents. Due 
to the presence of unknown lattice renormalisations, two conditions are required. The 
most suitable choice in the case of the vector current is 


iv = arctan (ytAnjvyX^senLVAj 


^(I/+P,-)=0. (18) 
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Indeed, for asymptotic times, the pion state dominates the matrix element^ and the 
condition reads 

(0|v;-^|vr-) = 0. (19) 

In case of the axialvector current we choose the condition^ 


= 0 


or asymptotically 

(0 |i+|p-) = 0. 

In terms of (|T^. d Eqs. (|TH|). (PH) admit the solution 


tanwy 




( 20 ) 


( 21 ) 


( 22 ) 


tan OJA 




(23) 


Eqs. d, (El), © and (d allow the numerical determination of oj and of the ratio 
ZaIZy- 

It is obvious that the definition of the twist angle in the lattice theory is subject 
to 0(a) ambiguities. Different choices of the parity-restoration conditions, including 
also the form of the lattice currents, result in different definitions of the twist angle 
differing by 0{a) terms. The situation of full twist corresponds to w = uoy = = 7r/2. 

Numerically it is most convenient to use ojv = 7r/2 as a criterion. The reason is that 
a safe determination of the twist angle is obtained in the asymptotic regime where the 
lightest particle dominates as intermediate state. This is the pseudoscalar state in the 
case of toy which, as one would expect, delivers a better signal than the vector meson 
in case of ua- Therefore we impose [snd 


or asymptotically 



p-) = o 


(Ol^xok ) = 0 


(24) 


(25) 


and denote with the corresponding value of for the given /r. 

^At small time-separations, due the 0{a) breaking of parity, intermediate states with “wrong” parity may 
still play a role. 

^In El El the use of the temporal component for the currents was proposed. This choice is however 
not optimal: a scalar state with positive parity dominates in this case the matrix element in the continuum 
limit, but at finite lattice spacing the 0{a) breaking of parity introduces contamination by pion intermediate 
states which eventually dominate for light quark masses. 





Another possible determination of u)v is obtained by replacing in the currents 
with their divergences. For simplicity, we consider the case of the conserved vector 
current which avoids the introduction of a renormalisation constant: 


cot Cbv = i 




tJ. x/i - y / 


m 


PCAC 


Py 




(26) 


Here in the last step laiEi the Ward identity for the conserved vector current 


d;vY, = 2i^iP+ 


(27) 


and the definition © of the “untwisted” PCAC quark have been used. If the 

local vector current defined in Eq. © is used for the determination of coy instead of 
the conserved one, in Eq. (El the introduction of the renormalisation constant Zy is 
required. In this case one has 


cot U!y = i 




Py 


= Zy- 


mt 




(28) 


where Zy is determined as explained in the next subsection. Using the definition m 
for LOy one arrives at the following relation involving this time the twist angle uj: 


cotcu = Za — - - 

P 


(29) 


Notice that the factor Zy cancels in this relation which is, therefore, independent of 
the choice for the vector current employed for the determination of the twist angle uj. 

One can simply show that the two determinations of ujy given by Eqs. El and 
(El coincide under the assumption that the ratio of the correlators is independent 
of the time separation; this is in particular true for asymptotic times where the pion 
dominates. 

To have an effective automatic 0{a) improvement, meaning without large O(a^) 
effects, the critical line (/iKcr(o 5 p), p) has to be fixed in such a way that the lattice 
definition of the untwisted quark mass (e.g. defined above) is free, on that line, 

from mass independent 0(a) errors. For a definition of the critical line where this 
condition is not necessarily satisfied, one has to make sure that p > aA^. 

The issue of the choice of the critical untwisted mass has been raised by the work of 
Aoki and Bar jSB] and by the numerical results obtained in ISl. This problem has been 
further analyzed in several aspects [HEEIISE]. In ElIinilHS] the theoretical framework 
is twisted mass chiral perturbation theory (tmChPT) |25j where the cutoff effects are 
included in the chiral lagrangian along the lines of El El The works ESHH agree on 
the fact that choosing the critical mass by imposing = 0 (or ujy = 7r/2) allows to 

have automatic 0{a) improvement down to quark masses that fulfill p ~ a^A^. In [SS] a 
Symanzik expansion was performed (in an approach different from that of refs. issnn. 
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cf. jin] for a discussion) confirming the results of [HHl • For a discussion of these 
issues in numerical studies within the quenched approximation see jUilSlilj and the 
review HH. 


3.2 Determination of Z 


V 


We adopt here the procedure well known in QCD which relies on the non-renormalisa- 
tion property of the conserved current 14^ |H8| . A possible determination of Zy in 
TMQCD is given by 

^(1) ^ (0|t^to,ok~) 

^ (0|fo+=o,oK-) ■ 

Note that in TMQCD the time component of the vector current couples the vacuum 
to the pseudoscalar particle: in the most interesting region near full twist this coupling 
is maximal. (Note that at a// = 0 the analogous procedure has to rely on the noisier 
matrix element with the vector particle or on three point functions.) Alternatively 
Zy can be determined without direct use of the conserved current by exploiting the 
(exact) Ward identity for the vector current. This implies 

—2in 


(0|W=o,c 


TT ) = 




-(0|FV=ok" 


Inserting the above relation in a second determination of Zy is obtained; 

(0|P^to|7r-) 


7 ( 2 ) _ 

Zy - 


(o|i4=o,ok' 


(31) 


(32) 


Zy^ and Zy^ (differing by 0(a) terms) are mass dependent renormalisations. We 


obtain a mass independent determination of Zy by extrapolating Zy' to full twist 
= 0)- In this situation the theory is 0(a) improved and the Zy^ deliver an 
estimate of Zy with O(a^) error (also including terms). 


3.3 Physical quantities 

The knowledge of the twist angle cv allows the derivation of physical quantities of 
interest in QCD for a generic choice of Let us consider the case of the 

quark mass and the pion decay constant. It is convenient inniiiai22i here to use the 
conserved vector current since it possesses already the right continuum normalisation. 
The physical PCAC quark mass can be obtained from the Ward identity for the 

physical axialvector current: 

(a;)i+P-) = 2am^'=^^(P+p-) . (33) 

We use Eq. 0 in order to eliminate in m for w / 0 

cot UJ + Cab (sfo w)“l (34) 
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and insert the result in the Ward identity (El using isospin invariance for As a 
result we obtain: 

^j^PCAC _ 

^ 2sinti; {Px Py) sintu ’ 

where in the last step we used once again the Ward identity ((23). Inserting Eq. El 
into the last expression in the above equation, we arrive at the following relation for 
the untwisted quark mass 



m^cAc _ coscj . (36) 

In the remainder we shall also make use of a definition of the untwisted quark mass 
which already incorporates the renormalisation factor of the axial current: 


mPCAc ^ ^pcAc ^ ^ . (37) 

Analogously, for the physical pion decay constant /jr we use 

af^r = (arn^)“^(0|A+^Q o|7r“) = -f(om^ sina;)“^(0|t7,lo^o|7r“) . (38) 


Also here the matrix element on the right hand side can be replaced by the matrix 
element of the pseudoscalar density as in m giving 


o/tt 


— 2a// 
sin w 


(0|P+=o|7r-) . 


(39) 


Let us note that here the normalisation of /jr corresponds to a phenomenological value 
~ 130 MeV. If the local vector current is used in El instead of the conserved one, a 
factor Zy is missing: 


a/i;7r =-f(am^ sin a;) ^(0|I/lo^ol^ )> fvn = Zy^U 


(40) 


3.4 Results 

In Fig. 121 the local determination of ujy and oja is shown as a function of the time 
separation for a specific simulation point at positive untwisted quark mass. The nu¬ 
merical values of the twist angles ujy, to a and u are reported in Table |11 Notice that 
the simulation point at /3 = 0.74 and k = 0.159 is almost at full twist. 

Figs. Eland El show the determinations of /x^cr by extrapolating and cot ojy 

to zero. The theoretical dependence of the twist angle upon the untwisted bare quark 
mass //k can be obtained [HI by starting from the equation E3 

cotw = —— + 0{a) (41) 

fJ-R 


11 






where fin and m^Fi are the renormalised twisted and untwisted quark masses in the 
continuum limit 


— O' (Mk l^Kcr) • (43) 

Observe that the relation m holds up to 0{a) terms because the right hand side of 
the relation corresponds to a different definition of the twist angle compared to the one 
given in Section |S3 The two definitions only coincide in the continuum limit. 

By using the first of Eqs. m one obtains for ujv m 

cot UJV = {Zovlj)~^ (^K - /^Kcr) + 0{a) (44) 

ZoV = ZgZAZp^Zy^ . (45) 

Note that the angular coefficient of the linear fit gives the finite combination of renor¬ 
malisation factors Zov- Using as an input the determination of Z^lZy in Eq. (1171) one 
can obtain from this the combination ZpjZg. 

We use Eq. (Hill for a linear fit to /U^cr and Zqv^ see Table [HI for the results. As 
expected from the discussion in Sec. EH the condition = 0 gives results very 

close to those from the parity-restoration condition cotcuy = 0. We conclude that 
the two methods are essentially equivalent also from the numerical point of view. A 
discrepancy is observed between the extrapolation from positive and negative quark 
masses for the simulation point {3 = 0.67: we interpret this as a residual effect of the 
first order phase transition at the given value of the lattice spacing. (Whether first 
order phase transition or “cross-over” can only be decided in a study of the infinite 
volume limit.) Observe also that the Zqv comes out different for the two different 
signs of the quark mass: this is due to the breaking of symmetry under reflection of 
the untwisted quark mass induced by 0{a) terms [1^1 • The numerical discrepancy 
shows that these 0{a) corrections are relevant. An 0(a)-improved estimate of Zqv 
is simply obtained by averaging the determinations for negative and positive quark 
masses, corresponding to a Wilson average for the quantity under study. An analogous 
observation can be done for other combinations of renormalisation constants (see the 
following). 

Table El reports the determination of the renormalisation constants of the vector 
and axialvector currents Zy and Za- The ratio ZaIZv comes from the analysis of the 
the twist angles, Eq. (HZl). Using the direct estimate of Zy by Eq. (PH) we can also 
determine Za- Observe that the full twist extrapolations of ZAjZy from the two quark 
mass signs present large discrepancies, which in this case cannot be attributed to 0{a) 
effects (these should disappear at full twist). A possible explanation of the discrepancy 
could reside in the relatively bad quality of the data in the negative mass region. The 
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discrepancies in Za and ZpjZs are a consequence of that for In the light of 

these considerations we rely on the determinations for positive quark masses. 

The full twist extrapolations of Zy are shown in Figs. [S] and |H1 the values from 
the two signs of the quark mass are rather close, compatible with each other within 
statistical uncertainty. For the case (3 = 0.74 the extrapolation is very short, see TableEl 
for the numerical values with comparison with one-loop perturbative estimates m 
Table |7| also includes the determinations of the ratio ZpjZg from Zov (see Eqs. (IHl) . 
(gSl)). This quantity is of particular interest for simulations m of the theory with an 
additional mass-split doublet describing the strange and charm quarks m Defining 
res as the mass-ratio rad mg-, the positivity of the fermionic measure in the strange- 
charm sector imposes 

Zp ^ T es 1 
Zs rcs +1 

The most stringent condition considering the experimental bounds mi for nis and nic 
is 

^ > 0.89 . (47) 

Zs 

Our results and the tadpole improved perturbative determinations for ZpjZs (for 
Nj = 2) seem to indicate that already at our values of (3 this condition is satisfied. 

The results for the physical PCAC quark mass and pion decay constant /tt obtained 
from Eqs. (|35|) and (|88|) are listed in Table |H1 In Eigs. Eland|H1 the pion decay constant 
is plotted as a function of the quark mass. The simulation points for negative quark 
masses are not taken into account in the present discussion. The figures also include 
the determination of /tt by the axialvector current a formula similar to Eq. 

applies in this case where, however, the factor l/sino; is replaced by l/coscu. In the 
interesting region near full twist this introduces large fluctuations in the estimate of 
/jr, as one can see from the figures. Moreover in the case of the axialvector current, 
the decay constant has not yet the right normalisation of the continuum: a Za factor 
is still missing. On the contrary, in the case of the conserved vector current /j^ has 
automatically the physical normalisation rnmniiiij. If we exclude the lightest point 
at /3 = 0.67, which is likely to be under the influence of residual metastabilities, /tt 
seems to be characterised by a linear dependence upon the quark mass. On the basis 
of this observation we try a simple linear extrapolation to the chiral limit = 0, 

see Table El for the numerical results. Of course, deviations from this linear behaviour 
could be present for lighter quark masses where chiral logarithms play a role. 

In order to check the scaling between the two (3 values we need to fix the lattice 
spacing. This can be accomplished by extrapolating the value of rg to = 0. Also 

in this case we obtain two different values for the two different signs of the untwisted 
quark mass, again due to 0{a) effects. As for Zqv we take the average of the two 
values, which delivers an 0(a)-improved estimate of rg in the chiral limit. The results 
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are reported in Table El We obtain for the lattice spacing (assuming rg = 0.5 fm): 
a(0.67) = 0.1757(41) fm, a(0.74) = 0.1326(70) fm. Denoting the zero quark mass limit 
of the pion decay constant by 


/o 


lim 

mPCAC=Q 


u 


(48) 


we obtain for /oro: /oro(0.67) = 0.333(10), /oro(0.74) = 0.274(20). These values are 
not far from the phenomenological value {foro)phen = 0.308. (The errors here are only 
statistical. Systematic errors of the chiral extrapolation are not included.) 


4 Fits to chiral perturbation theory 

Chiral perturbation theory (ChPT) is an expansion around the limit of massless quarks 
in QCD pp. It describes the dependency of physical quantities on the quark masses 
in terms of expansions in powers of quark masses, modified by logarithms. In nature, 
however, quark masses have fixed values. The question of how observables depend on 
them functionally is experimentally unaccessible. Lattice gauge theory, on the other 
hand, offers the possibility to vary quark masses. Therefore it represents the ideal 
field of application of chiral perturbation theory. On the one hand, chiral perturbation 
theory allows to extrapolate results from numerical simulations of QCD into the region 
of small physical values for the up- and down-quark masses. On the other hand, lattice 
QCD can provide values for the low-energy constants of chiral perturbation theory. 

In chiral perturbation theory the effects of the non-zero lattice spacing a can be 
taken into account in form of an expansion in powers of a[2l I45LI46LI47L I48| . For the 
case of the Wilson twisted-mass formulation of lattice QCD this has been worked out 
in next-to-leading order in 113 m muni. 

The major purpose of the present paragraph is to provide a set of formulas derived 
from lattice chiral perturbation theory that can be used to analyze physical quantities 
such as the pion mass, decay constants and amplitudes. The novelty here is that these 
quantities have to be described across or nearby a phase transition. 

The ChPT formulas are expected to be applicable at sufficiently small values of 
the lattice spacing and quark mass. It is thus far from obvious whether the data 
obtained with the DBW2 action in this work can be described by them, hence it is 
interesting to confront the simulation data at our quark masses and lattice spacings 
with these formulas. Let us emphasize that we consider this investigation mainly as a 
methodological study that does not aim to extract physical values of the low energy 
constants in the first place. 

Properly determined parameters of the ChPT formulas in the continuum limit are 
independent of the lattice action. The parameters describing the dependence on the 
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lattice spacing do, however, depend on it. Therefore, in an Appendix we also present 
ChPT fits of some simulation data obtained previously with the Wilson plaquette gauge 
action M- 

The quark masses in chiral perturbation theory always appear multiplied by 2Bq, 
where Bq is a low-energy constant. A connection to lattice regularisation can be 
established by considering the renormalised quark masses defined in Eqs. (|H, (ESI) 
and 

PCAC _ ^A PC AC 

^xR ~ Zp ^ 

A common renormalisation factor \/Zp in and can be absorbed into Bq. 

However, since the multiplicative renormalisation of and ^ differs by a factor 

Za, this has to be taken into account when fitting lattice data (see below). 

The lattice spacing enters chiral perturbation theory in the combination 

p = 2Woa , ( 50 ) 

where Wq is another low-energy constant. 

For the low-energy constants of lattice QCD in next to leading order gniiiHi with 
two quark flavours we use the notation 



T 54 — 2 L 4 -|- L 5 , Lge — 2Le + Lg , W 54 — 21^4 -|- H 5 , Wse — 21^6 -I- Wg , (51) 

W = ^{Ws6-2Lse), W' = ^{W^e-W86 + Ls6), W = . (52) 

Experience in untwisted lattice QCD shows m that lattice artefacts are consider¬ 
ably reduced when observables are considered as functions of the PCAC quark mass 
instead of the renormalised lattice quark mass. (A possible reason is that the PCAC 
quark mass reabsorbs leading order 0(a) effects.) Therefore, in our case, instead of 
using m^p as a variable, we re-expand the physical quantities in terms of the PCAC 
quark mass in the twisted basis . Including the relevant prefactor we define 

(53) 

For the purpose of fitting data at constant p it is convenient to define the combination 

X = 2Bo^{mll^^f + pl. (54) 

(The attentive reader is certainly realising that we use the symbols x for different quan¬ 
tities. Nevertheless, both the notation for the fermion field of twisted-mass fermions 
and the mass parameters in ChPT are standard in the literature and we do not want 
to change neither of them in this paper.) Then, for the charged pion masses, chiral 
perturbation theory at next-to-leading order including lattice terms of order a gives 

ml^ = X+ 32^2 j 72 ^ -k ^{(-L54 -F 2L8q)x^ + 2{W - W)p XpcacI • (55) 
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Similarly for the pion decay constant and the one-pion matrix element of the pseu¬ 
doscalar density: 


Fo ^ 167r2F2 



XpCAC 

X 


(56) 



327r2Fo2 


1 


X In ^ + 4L86)x + (4Vh - 2W)p 



In the ChPT formulas the pion decay constant at zero quark mass (Fq) appears. In the 
conventional normalisation its phenomenological value is Fq ~ 86MeV. This is related 
to /o ~ 122 MeV used in the previous section by Fq = /o/\/2. Similarly, F^- and Gjr 
denote the pion decay constant and the one-pion matrix element of the pseudoscalar 
density, respectively, in this normalisation convention. 

The renormalisation scale A appearing in the one-loop contributions is taken to be 
A = AttFq as usual. Taking into acconnt the renormalisation factors, when nsing these 
expressions for fitting the lattice data, one writes 




(58) 


X = 


where B = BqZaIZp. 

4.1 Fit procedure 

For fitting the data as a fnnction of Eqs. (ESI, El and m are going to be used. 

The data for m.^-, F^r and as well as that for are afflicted with numerical 

errors. Therefore, a fit procedure has to be used which takes into account errors in 
both coordinates. The method with effective variances m treats the coordinates on 
unequal footings bnt is numerically not so convenient. We have decided to use the 
more appropriate method of generalised least-squares fits E2- 

Consider a data set containing N “measured” values for each of the D variables. 
They are collected in the vector y = (y i ,..., yN), where each element yi is itself a 
column vector with D elements yi = {yi, j = 1, ..., F. The true values for each data 
point, which have to be estimated together with the parameters, will be collected in 
the same way in a vector x = (xi,..., xn) with entries Xi = {xij} , j = 1,..., D. Now 
the set of measured data points {yi,j} represents a single realization of an experiment 
which occurs with a probability given by a joint distribution called “likelihood”. The 
likelihood is specified by a multivariate normal distribution L with mean values given 
by the exact values x and a ND x ND covariance matrix a = , i,k = 





( 59 ) 
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The process of data analysis amounts to the constrained maximisation of this like¬ 
lihood through the estimation of the values of x based on the knowledge of y, where 
the constraints enter through the fit-functions. Instead of maximising L it is more 
convenient to minimise its negative logarithm. The only non-constant term is given by 

L' = ^ (x - y) a-^ (x - y)'^ . (60) 

The fit-functions are given by a number F of model-functions Gi, which can be 
incorporated as, generally nonlinear, constraints on the relationship between the exact 
values collected in x. These functions also depend on a set of P parameters a = 
(ai ,... ,ap), whose values are to be determined. They can be written in the compact 
form G (x,a) = 0 with the P-dimensional column vector G = (Gi ,... ,Gf)- 

Maximisation of the likelihood L under the constraints G (x, a) = 0 is now equiv¬ 
alent to the unconstrained minimisation of C given by 

£ = i (x-y)(T"^ (x-y)'^-h AG, (61) 

where A is the P-dimensional row-vector of Lagrange multipliers. We implemented the 
minimisation of C using the Maple algorithm NLPSolve, which is based on routines 
provided by the Numerical Algorithms Group (NAG). 

In the present case the N different points of measurement correspond to different 
values of the hopping parameter k, which are completely independent of each other. 
Therefore we can assume the covariance matrix to be diagonal, 
where Ayjj- denotes the statistical error of yij. 

The errors of the model-parameters Oi are calculated using a Monte-Garlo approach. 
In K steps of an artificial Monte-Carlo procedure a new set of normally distributed 
values ^ is generated using the values of {vij} as means and 

aij as the variances. Now for every k an independent estimate for the parameters is 
calculated yielding in each step. Finally the errors Aai are given by the standard 
deviation of the set of > k = I,..., K . 

4.2 Results 

At P = 0.67 and at /3 = 0.74 results for Pj^, Gt^ and are available both for 

non-vanishing and for vanishing twist mass y. At /r = 0 only part of the data, namely 
for mo — rricr > 0, is reliable and is being used. 

By using the results in Table 1^1 for the values of r^/a extrapolated to the chiral 
limit, we express all quantities in units of MeV. For the value of the Sommer scale we 
assume tq = 0.5 fm = (394.6.. MeV)“^. This allows us to compare and to combine the 
results from the different values of p. 
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It is important to observe that the lattice spacing a(/3) is obtained from extrap¬ 
olation of To/a to the chiral limit. In presence of both positive and negative masses 
we take the average. This is a strong constraint on the fits, since the data have to 
reproduce the scaling behaviour dictated by vq. If the purpose is the determination of 
the low-energy constants, matching ratios like with ChPT would be preferable. 

However, in this exploratory study, we find it interesting to check that the scaling 
behaviour of different quantities is indeed consistent. 

We made combined fits of the three quantities as functions of for both values 

of f3, including lattice terms of order a. For the pion masses the expressions for the 
O(a^) lattice terms are known, but they cannot be fitted meaningfully. The value of 
Za, entering the fit functions, has been taken as input from the Monte Carlo data. As 
it varies with (5, we denote the corresponding values Za{P). The fits include data both 
for non-zero and zero twisted mass /r. 

The low-energy constants resulting from the fits are shown in Table 11171 In the first 
case data points with both positive and negative values of are fitted, whereas 

in the second case only those with > 0. (This latter choice corresponds to 

the procedure in Section El where also only points with > 0 have been taken 

into account in the chiral extrapolation of f^-) We also made single fits for the three 
quantities, which are not displayed here. As they are each based on less data, their 
results are less valuable, but consistent with the combined fits. 

In addition to the single-/? fits we also made a global fit including the data from both 
values of /?. The results are also contained in Table Uni The fits at the two single values 
of (3 and the global fit are roughly consistent with each other. The differences in the 
numbers for the low-energy constants give an indication of the size of the uncertainties. 

Instead of using Za as input from the numerical calculations, it can alternatively 
be left as an additional fit parameter. The corresponding fit results are shown in the 
right hand side of the table. The fitted Za is in rough agreement with its Monte Carlo 
estimate. Also, the low-energy coefficients are consistent with the ones from the other 
fits. 

In addition to the combinations of Gasser-Leutwyler coefficients L^., the table in¬ 
cludes the values of the invariant scale parameters P] 

As = AttFq exp(1287r^(L54 - 21.86)), 

A 4 = dvrFo exp(327r^L54). (62) 

The results for A 3 , A 4 are close to phenomenological estimates (see the discus¬ 
sion). The W parameters have large errors but their magnitude is reasonable, as Wq 
is expected to be of order Ag^^ and the other VF’s of the same order as the L’s. 

The fit curves for F^^ and G-,^ together with data points at /? = 0.67 and (3 = 0.74 
are shown in Figs. muniiii] and nn In order to display the size of the leading order 
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contribution and the corrections, the figures contain additional curves representing the 
fit functions with some of the low-energy constants being put to zero. 

We have also investigated as a function of rriQ. It can be fitted with the 

corresponding formula from chiral perturbation, which involves W, W' and W but no 
L-coefficients, but the resulting coefficients are unreliable owing to large errors. 

In this section we stick to the definition of the untwisted bare PC AC quark mass 
in Eq. (jH)). As it is shown in Appendix lA.II the agreement with the ChPT 
formulas can be improved by taking of Eq. dsn as the quark mass variable, 

instead. In addition, ChPT fits to some previously obtained simulation data by the 
Wilson plaquette gauge action are also presented there. 

5 Discussion 

We compared in this paper the numerical simulation results with two flavours of 
twisted-mass Wilson quarks and DBW2 gauge action at two values of the lattice spacing 
corresponding to /3 = 0.67 and (5 = 0.74. The lattices were 12^ • 24 and 16^ • 32, respec¬ 
tively. The lattice spacing was defined by the value of the Sommer scale parameter vq 
extrapolated to zero quark mass and assuming tq = 0.5 fm. The /3-values were cho¬ 
sen in such a way that the lattice extensions were approximately equal: L ~ 2.11 fm 
and L ~ 2.12 fm, respectively. Also the bare twisted masses scaled approximately: 
ro/i ~ 0.0285 and ro/i ^ 0.0283, respectively. 

The comparison of the observed quantities at two /3-values allows for a first look 
at discretisation errors. The outcome of these tests is reasonable, having in mind 
the coarse lattice spacings: a ~ 0.176 fm on the 12^ • 24 and a ~ 0.133 fm on 16^ • 
32. Eor instance, the results for the pseudoscalar decay constant at zero quark mass 
are /j^ro = 0.330(10) and f-j^ro = 0.274(20) at /3 = 0.67 and (3 = 0.74, respectively. 
These values also come close to the phenomenological value {f-wrQ)phen = 0.308 [HH]- 
The situation is somewhat worse for the pseudoscalar-vector mass ratio, as Figure ^ 
indicates. There are some noticeable scale breaking effects, especially for pseudoscalar 
masses near ttItj = Of course, one has to have in mind that the p-meson mass in 
most of the points is quite close to the cutoff. 

The prerequisite for the extraction of quantities as, for instance, /tt is the know¬ 
ledge of the multiplicative renormalisation Z-factors for the currents. For obtaining 
the Z-factors one can exploit the twist-angle dependence in the plane of untwisted 
and twisted quark masses. As we have shown in Section |S1 this is a rather powerful 
method for obtaining “finite” (according to perturbation theory) Z-factor combinations 
as Zy, Za and ZpjZs- Remarkably consistent results could be obtained even with our 
exploratory simulation data, without a dedicated choice of simulation points for this 
purpose. 
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We have also attempted to describe our numerical simulation data by a set of 
formulas derived from lattice chiral perturbation theory. Although the values of the 
lattice spacing and the quark mass are rather large in the simulations, it turned out 
that these formulas describe the behaviour of many physical quantities - even across 
the phase transition - surprisingly well, at least on a qualitative level. However, at the 
quantitative level our presently available data do neither allow to make a quantitative 
extraction of the values of the ChPT parameters nor can we answer the question 
whether the lattice artifacts are well described by the lattice extension of ChPT. The 
achieved qualitatively correct ChPT fits of our simulation data makes us very optimistic 
that with new data we are working on at present - at smaller lattice spacings and small 
quarks masses - these questions will be answered. To achieve this the experience with 
the fits in this paper will be very helpful. 

In Section 0] and Appendix IA.1I we used the NLO expressions of ChPT including 
terms describing 0{a) lattice artefacts. In general, metastable points near the first 
order phase transition can be and have been included in the fits. (Note that the fits in 
Ref. [S3 a-i'e also based on metastable points, as it has been discovered later.) Several 
setups were tried and were shown to give satisfactory and consistent fits. Nevertheless, 
there are probably some higher order effects (higher orders both in the quark mass and 
in lattice spacing) which are non-negligible in our parameter range. In addition, for the 
multi-parameter fits our data are not precise enough and the data points are too few 
and not optimally distributed in the parameter space. (In a dedicated investigation 
the inclusion of partially quenched data points could be very helpful.) Qualitatively 
speaking, the ChPT fits presented here support the choice of the PCAC quark mass as 
the preferred quark mass variable and show that the 0{a) effects are not overwhelming 
because a fit without them is most of the time possible. Both these findings agree with 
those of Ref. [201 • 

The ChPT fits are also helpful in estimating the minimal pion mass at a given 
lattice spacing. For instance, the results at /3 = 0.74 (a = 0.1326(70) fm) indicate that 
for fixed = 0.0075 we are above the end point of the first order phase transition line 
(see e.g. the smooth behaviour near in Figure jlj). The minimal value of the pion 
mass in Figures inil and IT^ is about m™"(a/r = 0.0075) ~ 280MeV. This is an upper 
bound for the absolute minimum at (3 = 0.74. 

According to Table El the fits of the data with DBW2 gauge action suggest the 
following qualitative estimates for the values of the relevant ChPT parameters: 

2.9 GeV < R < 3.5 GeV 
70 MeV < Fq < 85 MeV 
4.0 < A 3 /F 0 ^ 8.0 

16.0 < A 4 /F 0 < 19.0 (63) 
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As Table im shows, the fits with the plaquette gauge action are roughly consistent with 
these values. The estimates of A 3^4 are close to previous estimates in EO]: A 3 /F 0 « 8 , 
A4/F0 21. 

The values of the IT-parameters describing 0(a) effects are not well determined 
and are in most cases consistent with zero in our fits, if (or is taken 

as the independent variable. Note that if one considers the relation of arn^^'^ and 
arriQ then W and W are quite visible. An example is Figure 2 in our previous pro¬ 
ceedings contribution EH] where W gives the difference of the slope between positive 
and negative masses (W' comes out to be small). 

If the data at the two /3-values are fitted separately, as Table 11171 shows, there is a 
remarkably good agreement of the corresponding parameter values. This agrees with 
expectations since the inclusion of 0{a) terms in the formulas reduces the discretisation 
errors in the physical parameters. The consistency of the ChPT fits is supported by 
the agreement of the pion decay constant at zero quark mass Fq with the value directly 
extracted from the data in Section |3| /o(/3 = Q.7A)/y/2 ~ 76MeV. The estimates of 
the universal low energy scales A 3^4 are within the bounds of their phenomenological 
values given in [SS]: A 3 = 0.6 (-1-1.4,—0.4) GeV, A 4 = 1.2 (-1-0.7,—0.4) GeV that is 
2.3 < A 3 /F 0 < 23.3, 9.3 < A 4 /F 0 < 22.1. 
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A Appendix 

A.l Comparison with the fits to plaquette-action data 

It is interesting to compare the results obtained from the DBW2 gauge action with 
those presented in m resulting from the plaquette gauge action. As shown in the 
previous sections, Chiral Perturbation Theory for Wilson lattice fermions (WChPT) 
offers a natural framework to perform such comparison. In fact, if NLO WChPT is 
applicable, the parameters Bq, Fq and L, entering equations (EH), (EEl) and (EZl) should 
already take their physical (continuum) values: lattice artefacts are expected to be 
taken into account by the W parameters. The latter depend, in general, on the lattice 
action. 

We remark that, having expressed all quantities (TV, Gtt and m.,^) as functions of 
rn^cAc^ the parameter W' (see 13 CS]) disappears, and the pion mass can apparently 
go to zero when —> 0 and /i ^ 0. However, one should keep in mind that not 

all values of are accessible with stable simulation points. This parametrization 

allows to include in the ChPT fit also metastable points, where both and are 

lower than it would be possible in a stable minimum of the effective potential. Since 
this is an interesting check, we exploit this possibility and we include also metastable 
points (from M) in the fit. 

Given the larger amount of data points, we use a different fit procedure from the 
one described in Section sm The is defined as for the effective variances method 
EDISaiMl, but minimized through the Matlab implementation of the Nelder-Mead 
Simplex Method. The variables a and are taken as independent variables, and 

Ft^, G-jt and mjr as dependent ones. 

Besides using a fitting procedure different from the one in the previous Sections 14. II 
Ea our fits to the plaquette gauge action data are restricted to data points with non¬ 
zero twisted mass {a/i > 0) only. We also tried to use different independent variables 
instead of , which correspond to different possible definitions of the untwisted 

component of the quark mass. It turned out that the fit quality is improving if one 
considers defined in (IH7I) . The difference implied by these changes compared 

to the analysis in Sections 14.114.21 - i.e. different fitting procedure, restricting the fit to 
a/i > 0 and using - is illustrated by Figures El El and ESI which have to be 

compared to Figures 1^ [TUI (TUI and El respectively. 

A consequence of considering instead of is that enters only 

indirectly - through the determination of a; - therefore we do not need to fit them. As 
said before, the Zp is included in the B factor. However, when comparing different 
lattice spacings and different actions, we must allow a /? dependent Zp. In practice we 
choose a reference {3 (corresponding to the smallest a which appears in the fit) and we 
fit a correction to Zp for each different a. These are not given in the table, but they 


25 




are always between 0.95 and 1.35. 

We summarize our results for the plaquette gauge action data in Table ^2 No 
statistical errors are quoted, since the systematic errors dominate, as the comparison 
of the results from the different fit setups shows. We perform fits including all data 
(top part of Table [TT|) or only data at positive mass (bottom part of Table fTTI) . In this 
second case the W parameters are set to zero. 

In Figures dlllZI and m the fits of the plaquette gauge action data are presented. 
Similarly to the DBW2 fits, the W's are very unstable, depending on the chosen subset 
of data, and in general they are consistent with zero within errors. The physical 
combinations L 54 and Lgg are consistent with the values obtained by the DBW2 gauge 
action. 

We also performed fits of all data and imposing W = W = Q. The values of the 
physical quantities are still reasonable, however the curves fit the data worse. We have 
also attempted fits where all the NLO parameters are set to zero {Li = IT = IF = 0), 
or where only lattice artefacts are included {Li = 0). Both these assumptions result in 
very poor fits, essentially because they cannot reproduce the curvature in Ft^ and G^r- 
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Tables 


Table 1: Run parameters: the gauge coupling (f3), the twisted mass in 
lattice units (apt,), the hopping parameter (n) and the lattice size. The 
last column shows the number of gauge configurations used in the data 
analysis. 


run 

P 

apt 

n 

L3 

X 

T 

.^conf 

(a) 

0.67 

0 

0.1650 

12^ 

X 

24 

4514 

(b) 

0.67 

0 

0.1655 

12^ 

X 

24 

2590 

(c) 

0.67 

0 

0.1660 

12^ 

X 

24 

2589 

id) 

0.67 

0 

0.1665 

12^ 

X 

24 

1721 

W) 

0.67 

0.01 

0.1650 

12^ 

X 

24 

600 

(&') 

0.67 

0.01 

0.1655 

12^ 

X 

24 

620 

(cO 

0.67 

0.01 

0.1660 

12^ 

X 

24 

509 

id') 

0.67 

0.01 

0.1665 

12^ 

X 

24 

570 

ie') 

0.67 

0.01 

0.1670 

12^ 

X 

24 

584 

if) 

0.67 

0.01 

0.1675 

12^ 

X 

24 

499 

if) 

0.67 

0.01 

0.1680 

12^ 

X 

24 

606 

iA) 

0.74 

0 

0.1580 

16^ 

X 

32 

1319 

iB) 

0.74 

0 

0.1585 

16^ 

X 

32 

419 

iA') 

0.74 

0.0075 

0.1580 

16^ 

X 

32 

430 

iB') 

0.74 

0.0075 

0.1585 

16^ 

X 

32 

296 

iC) 

0.74 

0.0075 

0.1590 

16^ 

X 

32 

353 

iD') 

0.74 

0.0075 

0.1595 

16^ 

X 

32 

352 
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Table 2: The results for the scale parameter (ro/a), the pseudoscalar 
(“pion”) mass am.,, and the vector-meson (“p-meson”) mass amp. 


run 

rp/a 


amp 

m^/mp 

rom^ 

irpm^)^ 

(a) 

2.305(36) 

0.4468(30) 

0.7025(44) 

0.6359(51) 

1.030(19) 

1.061(38) 

ib) 

2.391(56) 

0.4085(55) 

0.7007(79) 

0.5831(66) 

0.977(23) 

0.954(44) 

(c) 

2.351(27) 

0.3619(27) 

0.629(10) 

0.5747(84) 

0.850(11) 

0.724(19) 

{d) 

2.652(38) 

0.235(12) 

0.595(22) 

0.396(18) 

0.623(30) 

0.389(37) 

{a') 

2.347(26) 

0.4540(24) 

0.7026(46) 

0.6461(47) 

1.065(12) 

1.135(25) 

ip') 

2.415(24) 

0.3981(40) 

0.6808(66) 

0.5847(61) 

0.9618(25) 

0.925(18) 

ic') 

2.503(29) 

0.3449(40) 

0.662(11) 

0.520(10) 

0.863(11) 

0.745(18) 

id') 

2.867(29) 

0.2793(26) 

0.654(45) 

0.426(30) 

0.801(16) 

0.641(26) 

ie') 

3.127(31) 

0.2937(32) 

0.807(64) 

0.363(29) 

0.918(14) 

0.844(25) 

if) 

3.279(36) 

0.3706(50) 

0.913(72) 

0.403(33) 

1.215(23) 

1.477(57) 

if) 

3.261(31) 

0.4514(84) 

1.013(82) 

0.444(36) 

1.472(30) 

2.168(88) 

iA) 

3.563(33) 

0.3038(15) 

0.5256(37) 

0.5780(41) 

1.082(11) 

1.172(23) 

iB) 

3.741(90) 

0.2250(29) 

0.491(14) 

0.457(13) 

0.843(22) 

0.711(36) 

iA') 

3.467(51) 

0.3107(24) 

0.5354(71) 

0.5803(78) 

1.077(17) 

1.161(36) 

iB') 

3.78(10) 

0.2429(36) 

0.537(21) 

0.451(18) 

0.920(25) 

0.846(46) 

iC) 

3.87(10) 

0.1954(22) 

0.57(14) 

0.337(79) 

0.756(31) 

0.572(48) 

iD') 

4.148(65) 

0.2620(38) 

0.639(73) 

0.409(48) 

1.086(24) 

1.181(52) 






































































































































Table 3: The results for the PCAC quark mass (arn^^'^) and pseudoscalar 
(“pion”) decay constant (af^T,). 


run 

O 

> 

O 



^o/xTT 

(a) 

0.03884(22) 

0.0895(14) 

0.18567(90) 

0.4279(62) 

ib) 

0.03224(71) 

0.0771(18) 

0.1798(17) 

0.4301(98) 

(c) 

0.02247(80) 

0.0528(20) 

0.1553(27) 

0.3653(75) 

{d) 

0.00972(43) 

0.0258(11) 

0.1369(65) 

0.363(18) 

W) 

0.03801(63) 

0.0892(16) 

0.05774(88) 

0.1355(25) 

ip') 

0.02791(65) 

0.0674(16) 

0.0520(11) 

0.1257(28) 

ic') 

0.01846(99) 

0.0462(22) 

0.0442(20) 

0.1107(44) 

id') 

0.00505(82) 

0.0145(22) 

0.0174(26) 

0.0499(75) 

ie') 

-0.0109(2) 

-0.0341(37) 

-0.0354(37) 

-0.110(12) 

if) 

-0.0252(18) 

-0.0829(62) 

-0.0562(44) 

-0.184(15) 

if) 

-0.0409(17) 

-0.1336(56) 

-0.0683(30) 

-0.2229(98) 

iA) 

0.02313(23) 

0.0824(10) 

0.1243(12) 

0.4429(58) 

iB) 

0.01251(43) 

0.0469(24) 

0.1124(37) 

0.420(22) 

iA') 

0.02247(33) 

0.0779(16) 

0.03645(60) 

0.1264(28) 

iB') 

0.01093(49) 

0.0414(21) 

0.0266(12) 

0.1007(46) 

iC) 

-0.00120(18) 

-0.0046(29) 

-0.0043(18) 

-0.016(11) 

iD') 

-0.01635(66) 

-0.06783(29) 

-0.0361(16) 

-0.1500(71) 
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Table 4: The twist angles oj, ujy and oja, o-s defined in Eqs. and 

m. determined by Eqs. dm, (Hf and m. 



afji 

K 

Uv/tT 

UJa/t^ 

uj/n 

0.67 

1.0 

■ 10-2 

0.1650 

0.1352(13) 

0.0564(17) 

0.0883(13) 

0.67 

1.0 

■ 10-2 

0.1655 

0.1772(29) 

0.0771(27) 

0.1190(25) 

0.67 

1.0 

■ 10-2 

0.1660 

0.2412(62) 

0.1069(41) 

0.1661(54) 

0.67 

1.0 

■ 10-2 

0.1665 

0.411(12) 

0.229(17) 

0.334(17) 

0.67 

1.0 

■ 10-2 

0.1670 

0.678(12) 

0.622(16) 

0.647(11) 

0.67 

1.0 

■ 10-2 

0.1675 

0.8053(86) 

0.826(13) 

0.8137(80) 

0.67 

1.0 

■ 10-2 

0.1680 

0.8709(43) 

0.843(23) 

0.857(11) 

0.74 

7.5 

■ 10-3 

0.1580 

0.1542(26) 

0.0722(38) 

0.1076(31) 

0.74 

7.5 

■ 10-3 

0.1585 

0.2613(66) 

0.1393(96) 

0.1963(77) 

0.74 

7.5 

■ 10-3 

0.1590 

0.532(12) 

0.582(37) 

0.5544(92) 

0.74 

7.5 

■ 10-3 

0.1597 

0.7966(49) 

0.790(15) 

0.794(12) 


Table 5: Determination of la^cr by requiring eo = tt/ 2, Hncrii^v), or 
m^cAc _ . The plus and minus signs indicate extrapolations 

from positive or negative untwisted quark masses avg denotes the 

average. 



ap, 

sign 

f^Kcr ) 


ZoV 

0.67 

1.0 

■ 10-2 

+ 

2.99800(9) 

2.99839(12) 

1.438(33) 

0.67 

1.0 

■ 10-2 

— 

3.00059(13) 

3.00043(17) 

1.065(61) 

0.67 

1.0 

■ 10-2 

avg 

2.99930(11) 

2.99941(15) 

1.251(47) 

0.74 

7.5 

■ 10-3 

+ 

3.145528(52) 

3.145645(22) 

1.328(36) 

0.74 

7.5 

■ 10-3 

— 

3.145441(52) 

3.145435(21) 

1.055(49) 

0.74 

7.5 

■ 10-3 

avg 

3.145484(52) 

3.145540(22) 

1.191(42) 
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Table 6: Renormalisation constants of the vector and axialvector currents. 
The ratio ZaIZv is determined from the analysis of the twist angles, cf. 
Eg. (Egi; two different determinations of the vector current Zy are re¬ 
ported: Zy^ from Eg. /T^) and Zy^ from Eg. Wi the renormalisation 
constant of the axialvector current is derived by combining the results for 
Za/Zv and Zy\ 


/3 

ajji 

K 

Za/Zv 



Za 

0.67 

1.0 

■ 10-2 

0.1650 

1.589(26) 

0.5910(13) 

0.5810(16) 

0.939(15) 

0.67 

1.0 

■ 10-2 

0.1655 

1.587(28) 

0.5813(11) 

0.5761(25) 

0.923(16) 

0.67 

1.0 

■ 10-2 

0.1660 

1.649(28) 

0.5766(12) 

0.5708(38) 

0.951(16) 

0.67 

1.0 

■ 10-2 

0.1665 

1.979(68) 

0.5689(10) 

0.5657(39) 

1.126(39) 

0.67 

1.0 

■ 10-2 

0.1670 

0.815(58) 

0.5705(14) 

0.5666(46) 

0.465(33) 

0.67 

1.0 

■ 10-2 

0.1675 

1.087(47) 

0.5716(32) 

0.5688(38) 

0.623(27) 

0.67 

1.0 

■ 10-2 

0.1680 

0.894(78) 

0.5851(33) 

0.5754(43) 

0.518(46) 

0.74 

7.5 

■ 10-3 

0.1580 

1.508(35) 

0.6379(12) 

0.6315(32) 

0.963(22) 

0.74 

7.5 

■ 10-3 

0.1585 

1.515(59) 

0.6294(11) 

0.6294(38) 

0.953(37) 

0.74 

7.5 

■ 10-3 

0.1590 

1.65(45) 

0.62595(95) 

0.6241(38) 

1.04(28) 

0.74 

7.5 

■ 10-3 

0.1597 

0.972(73) 

0.6291(25) 

0.6242(40) 

0.612(46) 
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Table 7: Full twist extrapolations for Zy, Za and the ratio ZAjZy (see text 
for explanations) with comparison with 1-loop perturbative estimates (FT) 
and tadpole-improved perturbative estimates (TI)M- The ratio Zp/Zs is 
also reported, determined from Zoy (see Eqs. M, m)- 


p 

apL 

Sign 

Op. 

z 

Z(PT) 

Z(TI) 

0.67 

1.0 

■ 10-2 

+ 

V 

0.5650(11) 

0.6089 

0.6531 

0.67 

1.0 

■ 10-2 

- 

V 

0.5673(19) 

0.6089 

0.6531 

0.74 

7.5 

■ 10-3 

+ 

V 

0.6217(23) 

0.6459 

0.6892 

0.74 

7.5 

■ 10-3 

- 

V 

0.6257(10) 

0.6459 

0.6892 

0.67 

1.0 

■ 10-2 

+ 

A 

0.952(30) 

0.7219 

0.7176 

0.67 

1.0 

■ 10-2 

- 

A 

0.49(4) 

0.7219 

0.7176 

0.74 

7.5 

o 

1 

CO 

+ 

A 

0.944(74) 

0.7482 

0.7735 

0.74 

7.5 

O 

1 

CO 

- 

A 

0.612(46) 

0.7482 

0.7735 

0.67 

1.0 

■ 10-2 

+ 

A/V 

1.683(52) 

1.1130 

0.9696 

0.67 

1.0 

■ 10-2 

- 

A/V 

0.867(70) 

1.1130 

0.9696 

0.74 

7.5 

O 

1 

CO 

+ 

A/V 

1.52(12) 

1.1023 

0.9747 

0.74 

7.5 

■ 10-3 

- 

A/V 

0.972(73) 

1.1023 

0.9747 

0.67 

1.0 

■ 10-2 

+ 

P/S 

1.17(6) 

0.8157 

0.9407 

0.67 

1.0 

■ 10-2 

- 

P/S 

0.81(11) 

0.8157 

0.9407 

0.74 

7.5 

■ 10-3 

+ 

P/S 

1.14(12) 

0.8302 

0.9444 

0.74 

7.5 

■ 10-3 

- 

P/S 

0.92(10) 

0.8302 

0.9444 
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Table 8: Physical PC AC quark mass and pion decay constant 

af^, obtained from Eqs. (ESI and respectively. The last two columns 
show = cos(a;) and the unrenormalised pion decay constant 

calculated with the local current af^K, respectively. 


/? 

ap 

K 

ajyjPCAC 



^fvTT 

0.67 

1.0 

■ 10-2 

0.1650 

0.03652(53) 

0.1672(25) 

0.03511(54) 

0.2936(63) 

0.67 

1.0 

■ 10-2 

0.1655 

0.02739(55) 

0.1541(25) 

0.02549(59) 

0.2750(73) 

0.67 

1.0 

■ 10-2 

0.1660 

0.02006(59) 

0.1447(23) 

0.01739(69) 

0.2549(84) 

0.67 

1.0 

■ 10-2 

0.1665 

0.01154(11) 

0.1192(18) 

0.00575(71) 

0.2113(62) 

0.67 

1.0 

■ 10-2 

0.1670 

0.01117(38) 

0.1085(37) 

-0.00497(43) 

0.1932(80) 

0.67 

1.0 

■ 10-2 

0.1675 

0.01810(69) 

0.1203(44) 

-0.01508(82) 

0.219(13) 

0.67 

1.0 

■ 10-2 

0.1680 

0.0230(17) 

0.1146(95) 

-0.0207(18) 

0.202(14) 

0.74 

7.5 

■ 10-3 

0.1580 

0.02262(45) 

0.1170(25) 

0.02133(66) 

0.1833(57) 

0.74 

7.5 

■ 10-3 

0.1585 

0.01297(44) 

0.0999(26) 

0.01057(54) 

0.1625(83) 

0.74 

7.5 

■ 10-3 

0.1590 

0.007611(38) 

0.0874(15) 

-0.00129(22) 

0.1400(56) 

0.74 

7.5 

■ 10-3 

0.1595 

0.01245(61) 

0.0867(39) 

-0.00992(78) 

0.137(10) 


Table 9: Chiral extrapolation = 0) of the Sommer scale param¬ 

eter vq and pion decay constant f^,. (This latter is denoted by /q = 
lim^pcAc=o The scale independent combination f^r^ is also reported. 
Only data with positive twisted quark masses have been used for the extrap¬ 
olations, with the exception of the point at ap. = 0.0075 and k = 0.1590 
which is almost at full twist. 


/3 

ap 

ro/a 

a [fm] 

o /o 

foro 

0.67 

0.74 

1.0 ■ 10-2 

7.5 ■ 10-3 

2.845(66) 

3.77(20) 

0.1757(41) 

0.1326(70) 

0.1171(59) 

0.0726(25) 

0.333(10) 

0.274(20) 
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Table 10; Results of the ChPT fits with DBW2 gauge action. Upper part: 
fit with both positive and negative values of . Lower part: fit with 

only positive values of . 



input Za 

fitted Za 


fi = 0.67 

= 0.74 

both 

= 0.67 

= 0.74 

both 

^a(0.67) 

0.952(30) 

- 

0.952(30) 

0.8658(90) 

- 

0.852(14) 

^a(0.74) 

- 

0.944(74) 

0.944(74) 

- 

0.868(18) 

0.909(31) 

Fo [MeV] 

80.7(3.6) 

68.6(5.2) 

73.7(4.8) 

78.9(3.2) 

66.0(4.4) 

72.0(3.0) 

5(0.67) [GeV] 

3.20(13) 

- 

3.20(12) 

3.09(10) 

- 

3.063(94) 

5(0.74) [GeV] 

- 

3.31(30) 

3.16(38) 

- 

3.12(19) 

3.18(15) 

L 54 ■ 10 ^ 

0.98(26) 

0.96(26) 

1.17(28) 

0.50(15) 

0.80(23) 

0.74(12) 

Lge ■ 10 ^ 

0.78(13) 

0.81(11) 

0.94(14) 

0.554(84) 

0.76(10) 

0.709(61) 

Wo-W- 10-3 [MeV3] 

50(15) 

-21(16) 

18(17) 

35(12) 

-30(14) 

6 . 6 ( 8 . 0 ) 


89(19) 

21(38) 

64(29) 

62(14) 

-9(22) 

35(12) 

A 3 / 5o 

6 . 1 ( 2 . 8 ) 

5.5(2.3) 

5.1(2.5) 

5.9(1.7) 

5.0(2.0) 

5.3(1.1) 

A 4/50 

17.1(1.4) 

17.0(1.4) 

18.2(1.6) 

14.74(69) 

16.2(1.1) 

16.86(59) 

^min! d.O.f. 

12.8(3.5) 

12.3(4.9) 

13.1(7.2) 

9.2(1.6) 

11.6(2.4) 

9.4(1. 6 ) 

Za{0.67) 

0.952(30) 

- 

0.952(30) 

0 . 888 ( 10 ) 

- 

0.896(11) 

Za{0.74) 

- 

0.944(74) 

0.944(74) 

- 

0.910(18) 

0.880(23) 

Fo [MeV] 

80.3(3.4) 

91.2(5.4) 

83.9(4.4) 

79.3(3.4) 

89.9(4.2) 

82.2(2.6) 

5(0.67) [GeV] 

2.92(11) 

- 

2.95(11) 

2.85(10) 

- 

2.864(84) 

5(0.74) [GeV] 

- 

3.46(22) 

3.52(38) 

- 

3.39(15) 

3.39(11) 

T 54 ■ 10 ^ 

1.39(33) 

1.04(53) 

1.32(28) 

0.86(17) 

0.82(23) 

0.80(13) 

Lge ■ 10 ^ 

0.92(16) 

0.71(20) 

0.81(15) 

0.70(11) 

0.64(13) 

0.649(77) 

A 3 / 5o 

7.1(4.!) 

7.7(6.4) 

7.7(4.0) 

6.4(2.2) 

6.9(3.0) 

6 .7(1.7) 

A 4/50 

19.5(2.0) 

17.4(2.9) 

18.6(1.6) 

16.47(88) 

16.3(1.2) 

16.19(69) 

d^min! d.O.f. 

10.1(4.9) 

2.7(6.5) 

5.8(7.9) 

7.0(2.0) 

2.6(2.4) 

4.5(1. 6 ) 
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Table 11: Results of the ChPT fits with plaquette gauge action. The 
columns correspond to different definitions of the currents for af^^ and 
. For the definitions see Section \d. lA Upper part: fit with both 
positive and negative values of . Lower part: fit with only positive 

values of . 




U & 



B [GeV] 

5.05 

5.04 

5.00 

4.90 

Fo [MeV] 

104.9 

104.2 

88.3 

86.6 

Ls6 ■ 10 ^ 

0.916 

0.950 

1.829 

1.943 

L 54 ■ 10 ^ 

1.637 

1.709 

2.850 

3.027 

ITo ■ IT ■ 10-3 [MeV3] 

31.5 

28.5 

2.9 

6.6 

ITo ■ IT ■ 10-3 [MeV3] 

43.2 

39.7 

-3.6 

-1.3 

A3/ Fq 

9.8 

9.9 

4.5 

4.2 

FjFo 

21.1 

21.6 

30.9 

32.7 

(Jfdev‘^la‘^)/d.o.f. 

2.08 

2.19 

4.25 

4.16 

B [GeV] 

5.05 

4.33 

4.43 

3.95 

Fo [MeV] 

98.5 

93.9 

90.5 

85.8 

Lse ■ 103 

0.892 

1.466 

1.135 

1.836 

F 54 ■ 103 

1.848 

2.705 

2.099 

3.155 

A3/ Fo 

13.6 

9.4 

10.1 

6.5 

A 4 /F 0 

22.5 

29.5 

24.4 

34.0 

(X] dev'^/a‘^)/d.o.f. 

1.36 

2.24 

1.26 

1.77 
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Figures 





Figure 1: The squared pion to p-meson mass ratio {rriT^/mp)'^ versus 
{rQUinY. Only simulation points with positive quark mass are considered. 
The physical point is shown by an asterisk. The straight line connecting 
the origin with it is the continuum expectation for small quark masses 
where both quantities are approximately proportional to the quark mass. 
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Figure 2: Determination of tan.ujv and tanujA as in Eqs. /fH) . (Hf for 
the point (a'). The lines represent the fitted values. 
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Rk= 

Figure 3: Determination of at 13 = 0.67, ajj, = 0.01 by parity- 
restoration and by extrapolating the untwisted PC AC quark mass 
to zero. 
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Figure 4: Determination of filler at 13 = 0.74, afi = 0.0075 by parity- 
restoration and by extrapolating the untwisted PC AC quark mass 
to zero. 


39 











0 0.01 0.02 0.03 0.04 


Figure 7: The pion decay constant af-j^ as a function of the PC AC quark 
mass at (3 = 0.67, au = 0.01. 









Figure 9: The charged pion masses squared as a function of at 

a/i = 0.01. The points represent the data at (3 = 0.67. The solid line 
displays the global fit with Za as input. The dashed and dotted lines show 
the fit with part of the L and W coefficients set to zero, in order to indicate 
the size of the NLO corrections. 
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Figure 10: The charged pion masses squared as a function of arn^^^^ at 
a/i = 0.0075. The points represent the data at 13 = 0.74. The solid line 
displays the global fit with Za as input. The dashed and dotted lines show 
the fit with part of the L and W coefficients set to zero, in order to indicate 
the size of the NLO corrections. 
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Figure 11; The pion decay constant aF^, as a function of at 13 = 

0.67, a/i = 0.01. The solid line displays the global fit with Za as input. The 
dashed and dotted lines show the fit with part of the L and W coefficients 
set to zero, in order to indicate the size of the NLO corrections. 


46 









Figure 12: The pseudoscalar matrix element a^G.,, as a function of 
at f3 = 0.67, ap, = 0.0075. The solid line displays the global fit with Za. as 
input. The dashed and dotted lines show the fit with part of the L and W 
coefficients set to zero, in order to indicate the size of the NLO corrections. 
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Figure 13: Fit of the charged pion mass squared from DBW2 data at non¬ 
zero a/i as described in Sec. lA.ll Upper (lower) curve belongs to j3 = 0.67 
(/3 = 0.74). 
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Figure 14: Fit of the pion decay constant aF^ from DBW2 data at non¬ 
zero an as described in Sec. lA.ll Upper (lower) curve belongs to j3 = 0.67 
(/3 = 0.74). 
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Figure 15: Fit of from DBW2 data at non-zero a/x as described in 
Sec. lA.ll Upper (lower) curve belongs to (3 = 0.67 {(3 = 0.74). 
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Figure 16; Fit of the charged pion mass squared from plaquette data 
at non-zero ap. Upper, intermediate and lower curves refer to (3 = 5.1, 
f3 = 5.2 and f3 = 5.3, respectively. 
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Figure 17: Fit of the pion decay constant from plaqnette data at non¬ 
zero a/i. Upper, intermediate and lower curves refer to f3 = 5.1, f3 = 5.^ 
and f3 = 5.3, respectively. 
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Figure 18: Fit of from plaquette data at non-zero a/i. Upper, 

intermediate and lower curves refer to f3 = 5.1, f3 = 5.2 and f3 = 5.3, 
respectively. 




